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Abstract 

We explore the role of cellular life cycles for viruses and host cells in an infection 
process. For this purpose, we derive a generalized version of the basic model of virus 
dynamics (Nowak, M.A., Bangham, C.R.M., 1996. Population dynamics of immune 
responses to persistent viruses. Science 272, 74-79) from a mesoscopic description. 
In its final form the model can be written as a set of Volterra integrodifferential 
equations. We consider the role of age-distributed delays for death times and the 
intracellular (eclipse) phase. These processes are implemented by means of prob- 
ability distribution functions. The basic reproductive ratio Rq of the infection is 
properly defined in terms of such distributions by using an analysis of the equilib- 
rium states and their stability. It is concluded that the introduction of distributed 
delays can strongly modify both the value of Rq and the predictions for the virus 
loads, so the effects on the infection dynamics are of major importance. We also 
show how the model presented here can be applied to some simple situations where 
direct comparison with experiments is possible. Specifically, phage-bacteria inter- 
actions are analysed. The dynamics of the eclipse phase for phages is characterized 
analytically, which allows us to compare the performance of three different fittings 
proposed before for the one-step growth curve. 

Key words: Virus dynamics. Cellular life cycle. Lytic cycle, Basic reproductive 
ratio 
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1 Introduction 

The interactions between viruses and cells in an infection process can be seen 
as an ecological system within the infected host. The mathematical description 
of these systems has attracted increasing interest in the last years (Wodarz, 
2006), especially concerning the characteristics of the immune response to a 
viral attack. A decade ago, Nowak and Bangham (1996) presented what has 
been called thereafter the Basic Model of Virus Dynamics (BMVD). This 
model has become quite popular among theorists and experimentalists (see 
Nowak and May (2000) and Perelson (2002) for some understanding reviews). 
The interplay between the BMVD and the effect of an immune response has 
proved useful to describe the dynamics of chronic HIV infections (Perelson, 
2002). Furthermore, it has provided interesting results regarding topics as the 
performance of drug therapies (Bonhoeffer et al., 1997; Wodarz and Nowak, 
1999), lymphocyte exhaustion (Wodarz et al., 1998), etc. 

The BMVD describes the time evolution of non-infected cells {X), infected 
cells (Y) and viruses {V) by the system of equations 



dX 

= X-SX - I3XV 
dt ' 

dY 

= j3XV -aY 
dt 

dV 

= kY - BXV - uV. 
dt ^ 



The infection process is governed by the parameter f3, which determines the 
rate of successful contacts between the target cells and the viruses. Mortality 
terms for the three species are considered with constant death rates 6, a and u, 
respectively. The parameter k measures the rate at which virions are released 



from a single infected cell. Finally, new target cells are produced by the host 
at a constant rate A. 

Despite the success achieved by the BMVD, it is clear that the model described 
in ([1]) is just a first approximation to the real underlying process. Probably 
the strongest simplification made in the model is that it assumes that the 
death rates are exponentially distributed (i.e., mortalities are considered as 
Markovian random processes) and therefore do not take into account accu- 
rately the details of the cellular life cycles. However, delays and structured 
life cycles are expected to play a very significant role in the dynamics of viral 
infections. For example, the infection process involves an intracellular phase of 
the virus, also known as the eclipse phase, which is not explicitly considered in 
([T]). For this reason, in the recent years some works have explored the effects 
of constant and distributed delays in the BMVD, also in the case where an 
immune response is considered. Herz et al., (1996) showed for the first time 
the importance of delays in order to explain the virus loads observed in HIV 
patients under drug treatment. This delayed model was later explored from a 
more formal point of view by Tam (1999). Similar ideas, with different expres- 
sions for the infection term, were considered by Culshaw and Ruan (2000), 
Fort and Mendez (2002) and Li and Wanbiao (1999). The effect of distributed 
delays was explored for different models of virus dynamics by Banks et al., 
(2003), Mittler et al., (1998) and Lloyd (2001). Finally, the role of a delayed 
immune response has been the subject of extensive research. Some examples 
are Buric et al., (2001), Canabarro et al., (2004), Wang et al., (2007) and the 
references there in, which focused on the chaotic patterns which can appear 
in these systems. 

The papers mentioned above have helped us to understand how delays can 



modify the cell-virus and virus- immune system dynamics. However, most of 
those works focused on the case where only one of the processes (usually the 
intracellular phase) is delayed. So, they do not considered the possibility of 
different delays for each process, whose combined contributions could modify 
the dynamical behavior of the system. 

On the other hand, the introduction of delays in the virus dynamics has been 
usually based on phenomenological (not always rigorous) arguments, without 
providing a justification of the delayed equations proposed. Only in Banks 
et al., (2003), Fort and Mendez (2002) and Wearing et al., (2005) a more 
formal discussion was provided. We stress that the implementation of delays 
into dynamical models is sometimes tricky, as memory effects can lead to the 
breakdown of hypothesis that are well established for Markovian processes. In 
fact, there is currently a very active research on this subject from the point 
of view of statistical mechanics (see, for example, AUegrini et al., 2003, Alle- 
grini et al., 2007, Rebenshtok and Barkai, 2007 and the references therein). 
According to these ideas, a rigorous mathematical approach is necessary to 
reach an accurate physical description of virus dynamics with delays. Here, 
we propose a system of Volterra integrodifferential equations which is a gen- 
eralization of the BMVD. This system of equations is derived from a meso- 
scopic approach where balance equations for each species {X, Y and V) are 
considered explicitly. Mesoscopic descriptions as that considered here (based 
on Continuous-Time Random Walk processes) have become quite usual tools 
for the description of physical and biological processes. At this stage, they 
have proved useful for the study of heat transport (Emmanuel and Berkowitz, 
2007), biological invasions (Mendez et al., unpublished), tumor cell growth 
(Fedotov and lomin, 2007), solute transport in porous media (Berkowitz et 
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al., 2000), earthquakes dynamics (Helmstetter and Sornette, 2002), financial 
markets (Masoliver et al., 2006) and many other. Here we will explore for the 
first time their application to the field of virus dynamics. 

Then, the aim of this paper is to use an integrodifferential approach to show 
how distributed delays can strongly influence the predictions from the BMVD. 
We find that the value of the basic reproductive ratio Rq and the values of the 
virus load can drastically change, in accordance with similar conclusions found 
in Lloyd (2001) from the analysis of the intracellular phase. Furthermore, 
the advantage of using such a general formalism as the one proposed here 
is that different situations of interest can be analyzed as particular cases of 
the model. According to this, we show how our model can be used to fit 
and characterize the one-step growth (osg) curve observed in phage-bacteria 
interactions. Three fittings proposed before by different authors are compared. 
We find that, albeit the three approaches fit reasonably well the osg curve, 
their predictions concerning the dynamics of the eclipse phase are slightly 
different. 

In the following, we show how a generalized version of the BMVD can be 
obtained using a mesoscopic description. In Section 2 we present our model, 
whose formal derivation is given in the Appendix for the sake of clarity. In 
Section 3 we explore the equilibrium states and their stability, which let us de- 
fine the basic reproductive ratio Rq. After that, we consider specific situations 
of special interest in virus dynamics. We consider the effects of distributed de- 
lays in the phase eclipse (Section 4) and in the mortalities for cells and viruses 
(Section 5). We also show how the model derived in Section 2 works in the case 
of phages-bacteria dynamics (Section 6), and we provide some examples using 
experimental data extracted from the literature. Finally, the main conclusions 
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obtained from our study are summarized in Section 7. 



2 The BMVD with distributed delays 

The model we consider here is depicted in Figure 1. It follows the same scheme 
as the BMVD but some of the random processes (those indicated by the dotted 
lines) are governed by their corresponding probability distribution functions 
(PDF). So that, (Pxit) represents the probability that a target cell X dies at 
age t, with equivalent definitions for v^y(t) and v^v'(t) for infected cells and 
viruses. Similarly, the function 0(t) determines the dynamics of the eclipse 
phase: a cell that becomes infected at time to can release 0(t) viruses at time 
to + t. 

The Volterra integrodifferential equations corresponding to the scheme in Fig- 
ure 1 read 



■-X-(3X{t)V{t)- I X{t-t')^x{t')^x{t-t',t)dt' 
Jo 

■■(3X(t)V(t)- I Y{t-t')^'Y{t')dt' 
Jo 

-(3X{t)V{t) + / (3X{t - t')V{t - t')0(t')$y(t')dt' 



dX{t) 
dt 

dm 

dt 
dVjt) 

dt i-'y^y y^y ■ j^ 

- f V{t - t')-^v{t')^v{t - t' ,t)dt' . (2) 

Jo 

The formal derivation of this model in terms of a mesoscopic description is 
provided in the Appendix. The functions "^x, ^y, ^y are defined by their 
Laplace transforms (we denote the Laplace transform of a function by the 
brackets [-js with the conjugate variable s) 

™^-[$xL ^"^^'-[^vl ^"^^'-[^vV ^^ 



where ^x{t) = J^ fxit')^^^' is the survival probabihty for the cells of age 
t. Analogous definitions hold for $y and ^y According to ([3]), the function 
\E'x(t) can be interpreted as the instantaneous death rate for a cell X of age t. 
Then, the term /q X (t — t')"^ x (t')^x (t — t' , t)dt' represents a generalized death 
term in which age-distributed death rates are considered, and where ^xit — 
t', t') is the probability that a particle X does not become infected during the 
time interval (t - t', t). Similarly, the term /q pX{t - t')V{t - t')(f){t')^Yit')dt' 
represents the release of new virions from those cells that became infected at 
time t — t', provided that these cells have survived up to time t. 

The system of equations ([2H3]) represents our generalization of the BMVD to 
the case with distributed delays. An important conclusion from ([2]) is that the 
density of infected cells Y does not appear in the equations for X{t) and V(t). 
It means that the formalism introduced here allows us to reduce the BMVD to 
a 2-species model. We do not need to consider explicitly the density Y{t)] the 
existence of the infected cells is implicitly considered by means of the function 
$y appearing in the equation for V{t). 

3 Equilibrium states and their stability 

The equilibrium states of the model ([2]) come from the analysis of the fixed 
points of the system at t ^ oo. There are two possible equilibrium states: the 
first one is the trivial, infection-free state, given by 

(x,„ye„v;,) = (Arx,o,o). (4) 

where we use Tj = /q°° $j(t)dt to denote the average lifetime of species i, with 
i = X,Y, V. The second state corresponds to the case of a successful infection 



defined by 



Xeq e-^^^'*<l>y(t)dt = — ^^^ ^^^ 

Jo Rq 

JO 



, e-Py^^^^^{t)dt = ^ (5) 

A 



where the equations fl27ll28p have been used, and we have defined 



Rq = (3Xtx 



■ roo 1 r /'CO 

/ e-^^"^*$y(t)dt / 0(t)<l>y(t)dt 

vo J \_Jo 



(6) 



As can be seen from ([5]), it is not possible to give exphcit expressions for the 
equihbrium densities. However, it can be proved that the infected state only 
has biological meaning {Yeq > and Veq > 0) if -Rq > 1- To see this, note that 
the condition i?o > 1 applied to the first equation of ([5]) implies Xeq < Xtx, 
which means that the equilibrium density in the infected state is lower than in 
the trivial state. Using that condition, it follows that the third equation in (^ 
has necessarily a positive solution for Veq. Hence, -Rq can be properly defined 
as the basic reproductive ratio, which is a key parameter in epidemiology and 
virus dynamics in order to predict the emergence of an infection (Anderson 
and May, 1991; Nowak and May, 2000). For Rq < 1 we have that every single 
virus generates statistically less than one new virus, so a permanent infection 
is not possible and the infected state does not exist. We also note that the case 
explored in the present paper, and so the expression ([6]), is more general than 
recent estimations for Rq where the possibility of a distributed intracellular 
period was also taken into account (Heffernan and Wahl, 2006). 

We will now explore the stability of the equilibrium states found. For this 
purpose, we will use the usual linear-stability analysis, so we introduce X{t) = 
Xeq + 6X{t) and V(t) = V^q + 6V(t). Inserting these definitions into ([2]) and 
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linearizing about the equilibrium state we obtain the following system for the 
perturbations 



^^^ = -m,SX{t) - (3X,,6V{t) - f 5X{t - t')^>x{t')dt' 
at JO 

+/3Xe„ 1 6V{t - t')^x{t')t'e-''^^'>''dt' 
Jo 

^^^ = ~l3X,,5V{t) - (3Ve,5X{t) + l3Ve, J^ 6X{t - t')<l){t')<^Y{t')dt' 

+(3Xeg f 5V{t - t')<p(t')^Y{t')dt' (7) 

JO 

- f 6V{t-t')mv{t')dt' + (3V,„ f 6X{t-t')mv{t')t'e-'^''^^''dt'. 
Jo Jo 

Since this system is now linear, we can propose for the perturbations expo- 
nential solutions of the form e'^* to get the characteristic equation 



= (/i + (3X,, + l^xl) (/i + f3X,, - (3X,g [<j)<l>yl + l^x] 

where we define [/]^ = J e~'^^f(t)dt in accordance with the notation used 
above for the Laplace transform. 

a) Infection-free equilibrium state 

First we analyze the stability of the trivial state corresponding to the absence 
of viruses. Introducing (jl]) into (jHD we obtain 

1 = PX,, [^vU,x., mrl ■ (9) 

From (Q, it is easy to find the necessary condition for the transition from 
stability to instability. In the BMVD it is known that the condition Rq ^ 1 
determines the stability of the infected- free state. From ([9]), it is possible to 
prove that, in general, this condition holds for any choice of the PDF's. The 
right hand side in that equation is a monotonically decreasing positive function 
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of /i and takes the value Rq a.t fi = 0. Then, if i?o > 1 both curves always 
intersect at a single point for a positive value of fi, which is nothing but the 
sufficient condition for the state to be unstable, independently of the PDF's 
considered. If i?o < 1 both curves always intersect at a single point but now 
for a negative value of /x. In this case the infection-free equilibrium state is 
linearly stable and infection dies out. 

b) Infected equilibrium state 

Using ([5]), the characteristic equation (jHj) for the infected state becomes ex- 
tremely complicated to treat, and it makes impossible to determine analyti- 
cally the stability of the infected state. However, we can still deduce the be- 
havior of this state by imposing some conditions to prevent the system from 
behaving unrealistically. First, we mention again that the infected state does 
not exist for _Ro < 1, so we only need to study the case -Rq > 1- Second, we can 
rewrite the first equation in ([5]), using Qj and the definition of the Laplace 
transform, as 

[0$y] - 1 

i^vU., = [^^^]^ ■ (10) 

Then, we conclude that there is only one possible positive solution for X^q, 
as the left hand side of this equation is a monotonically decreasing function 
of Xeq. From that, similar arguments can be applied to the third equation in 
([5]), so it follows that the solution for Veg is unique too. As a whole, we have 
that the infected state is always unique. This, together with the unstability 
of the non-infected state for i?o > 1, allows us to conclude that the infected 
state cannot be an unstable node or a saddle point, as it would imply that 
for some initial conditions the system would grow without control towards the 
state X -^ oo and/or V —>■ oo. This unbounded behavior is not possible in 
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our system. Then, the only possibihty is that the infected state is stable for 
Ro>l. 

The derivations presented in this Section show that the introduction of dis- 
tributed delays does not modify the stability conditions of the BMVD. Al- 
though our mesoscopic model ([2]) is much more general that the original ver- 
sion ([1]), we find that the condition Rq ^ 1 is always the one that determine 
the stability of the two possible equilibrium states. Note also that the condi- 
tion to have an infected state of coexistence between viruses and cells {Rq > 1) 
can be interpreted as a threshold value for the contact rate 

1 



P> 



Arx Jo°^e-^>^^x^^v{t)dt [/o°° 0(t)$y(t)dt] 
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4 The BMVD with a delayed eclipse phase 

We have presented a general model which takes into account distributed delays 
for the cellular death and the eclipse phase. However, the application of the 
general case requires knowing all the temporal distributions considered, which 
is not always possible at practice. Then, it can be useful to study some specific 
and simpler cases which have a special interest for application purposes. 

First, we consider the case where no age-distributed effects are introduced in 
the death process i.e. the probability of death is independent of the age of 
the cells. This corresponds to the situation used in the BMVD, which in our 
integrodifferential model is recovered by assuming ipx^ V^y? V^y ^s exponentially 
decaying functions (v^xl^) = Se^^^, V5y(t) = ae^"*, (fyit) = we""*). For the 
eclipse phase, we can assume that when a cell is infected, it takes a fixed 
constant time r until the first virion is released and after that, virions are 
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continuously released at a constant rate k. The delay r is the time necessary 
to inject the viral core into the cell and make its genetic machinery start the 
reproduction process. So that, the function 0(t) in our model will be taken as 
a step function (f){t) = kH{t — r), where H{) is the Heaviside function. 

This specific example has been studied by some authors before (Herz et al, 
1996; Tam, 1999; Culshaw and Ruan, 2000), so we can compare the predictions 
from our model with those previous approaches. Replacing the distribution 
functions (Pi{t), (pit) into the general model ([2]) we obtain 



^ = X-6X- (3XV 
dt ^ 

dV 

= pXV-aY 
dt ^ 

^ = t /3X{t - t')V{t - t')ke~''''dt' - (3XV - uV. (12) 

Q6 Jr 



In the equation for V{t), the expression j3X{t — t')V{t — t') represents those 
cells that became infected at time t — t' . So, the new virions appeared are 
equal to that expression multiplied by the rate k and by the probability e~"* 
that the infected cells have survived from time t — t' to t. The expression of 
Rq that one obtains for this case, from ([6]), is 



ii, = f Cle-" - l) . (13) 



Note that the system ( IT2|) is apparently different to the previous models pro- 
posed before for the analysis of a delayed eclipse phase (Herz et al., 1996; 
Tam, 1999; Culshaw and Ruan, 2000). In those works a delayed term f3X(t — 
T)V{t — t) was introduced ad hoc in the evolution equation for Y{t): 
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dX 

^ = X-6X- (3XV 

— = /3X(t - T)V{t - r)e-'^" - aY 

dV 

= kY-pXV-uV. (14) 



However, it is easy to see that the value of Rq for this model is exactly the 
expression flT3l) . and the equilibrium states coincide with those found from our 
model too. Actually, both models represent the same underlying process except 
for one subtle detail. In the model (fT4l) . the fraction of cells (3X{t ~T)V{t — T) 
are considered as infected cells only after the time delay r. But during the 
period from t — r to r these cells 'disappear', it is, they do not enter neither 
in the equation for Y nor in those for X or V. Instead, in our model the 
cells become Y cells at time t — t and they start releasing the new virions 
at time t, so our approach is phenomenologically more correct. Regarding the 
dynamics of both models, the only difference between (TT2|) and ( HM will be 
in the solution for Y{t): the value predicted by the model (fT4l) will be always 
below the real one, as some infected cells are not being counted. 



5 The effect of age-distributed times for cellular death 

Now we try to study a more realistic case according to the experimental data 
available in the literature. We will consider that the eclipse phase follows the 
same dynamics as that in Section 4. But the death times are now assumed to 
follow Gamma distributions, which are quite standard curves used for fitting 
experimental data to cellular death times (see for example the recent work by 
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Hawkins et. al. (2007)). Hence, in this case we will use 

0(t) = kHit - r) if^it) = ^^^l^a^^^^,^ (15) 

for i = X,Y, V, where r(-) denotes the gamma function and Oj and r* are the 
characteristic parameters of the Gamma distribution for mortality, with the 
average lifetime given by Tj = r*aj. 



Inserting these distributions into ([6]) the basic reproductive ratio Rq reads 

(16) 



o _ (l+^Ar^r^r-2 /., "^^ 



ay - 3 ( rV 



r- \n 



for ay integer. From flT6|) . it follows that the influence of distributed death ages 
could be important for the value of Rq and, as a result, it strongly modifies the 
value of the virus load at equilibrium. This effect is represented in Figure 2, 
which shows the numerical solution V{t) obtained from the model ([2]) for dif- 
ferent values of the parameter a (for simplicity we define a = ax = ay = ay). 
For a = 1 we recover the case where the death probabilities are exponentially 
distributed, it is, the prediction by the BMVD. In the three curves shown, 
the average lifetimes for the three species are kept the same. It allows us to 
compare properly the effects of the mortality distributions on the virus load 
dynamics. Two main differences are observed between the curves in Figure 
2. First, note that the virus loads decrease in time for t < 2; this is because 
we have used a delay r = 2 for the eclipse phase, so only after t = r the 
infected cells start to release the first virions, and then the virus load in- 
creases drastically. The minimum value observed at t = 2 is much lower in 
the case a = 1. This is because the BMVD assumes unrealistic high prob- 
abilities of death for the early stage of the infection, an effect which can be 
corrected by the Gamma-distributed mortalities used here. This point is of 
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great importance concerning the probabilities of a primary immune response 
to successfully clear the infection. Second, we also find important differences 
between the maximum virus loads reached at equilibrium; for the parameters 
used in Figure 2, the final virus load for a = 1 is approximately 10-fold higher 
than in the case a = 3. 

Therefore, we conclude that the BMVD underestimates the virus loads in the 
early stages of the infection and overestimates the peak of the virus load, if 
compared with the case of distributed mortalities considered here. In conse- 
quence, it turns out that we need to know with some detail the life cycle of 
viruses and cells to obtain an accurate picture of the infection dynamics. 

6 Application to phage-bacteria interactions 

The interaction between phages and bacteria can be described as two con- 
secutive steps: adsorption and reproduction (Mc Grath and Sinder, 2007). 
Adsorption consists in a collision between phage and bacteria resulting in a 
group, called infected bacteria, constituted by the bacteria and the phage at- 
tached to its membrane. The second step begins when the phage inoculates 
its genetic material into the host bacteria and begins to replicate it. From 
this time onwards the number of new viruses increases inside the bacteria, 
stopping when the bacteria bursts at the end of the latent period. Basically, 
the main difference between this situation and those explored in the previous 
Sections is that for phages the eclipse phase finishes with a lytic process that 
involves the death of the infected cell. In terms of the model presented here, 
this idea can be introduced simply by choosing the appropriate form for the 
function 0(t). 
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Here we deal with the reproduction process, which is known to produce a 
characteristic one-step growth curve V{t) for virulent phages. Let us consider 
that at t = the phage inoculates its genome and all the bacteria become 
infected instantaneously, with Y{t = 0) = Y{0). Then, we can define Jvit) = 
Y{0)(f){t) as the rate of viruses released at time t, following the same notation 
as in the Appendix (see Equation fl26p and the comments below). As all the 
cells are assumed to be already infected at t = 0, the infection process for 
t > can be obviated. We can thus take Qy = 1 in Equation (!26|) to obtain 

V{t) = V^(0)$y(t) + r r(O)0(t - t')<l>v{t')dt', (17) 

Jo 

which constitutes our theoretical model for the osg curve. If the osg is known 
from experiments, the function 0(t) can be determined by fitting that curve 
to some function and applying 

'dV r' ..„ ' 



^w-mU^/o''"-'''*^-"''* 



F(0) V dt _ 

^ ' \ / osg 

which comes directly from the solution of flT7|) . However, the result flTSj) can 
only be applied if we know the function \l/v', which is related to the mortality 
distribution ipy according to ([3]). At practice, the probabihty of death for the 
viruses is usually considered very small in the time scale of the experiments, so 
it can be neglected. In that case, "^v ~ and then we find that 0(t) becomes 
proportional to the derivative of the one-step growth (osg) curve 

For fitting the one-step growth V{t), some authors have considered before a 
piecewise function composed by three segments (You et al., 2002; Hadas et al., 
1997). Continuous functions have been proposed too, for example error func- 
tions (Rabinovitch et al., 1999) or logistic-like functions (Fort and Mendez, 
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2002; Alvarez et al., 2007). For these three cases one finds that the correspond- 
ing expressions for 0(t) are those shown in Table 1. We have written there the 
functions in terms of the parameters r, r and \4o- For the sake of completeness, 
we also show the relation between these parameters and the eclipse time, the 
rise rate and the burst size, which are commonly used in experimental works 
to characterize the osg curve (a proper definition of these is provided in Figure 
3). 

In Figure 4 we show with symbols the experimental results for one-step growth 
of phage T7 on E. coli BL21 grown at different rates (You et al., 2002), while 
the specific values obtained from the adjustment in each case are detailed in 
Table 2. The solid curves in Figure 4 represent the fitting of the experimental 
results to the logistic- like function, exhibiting a good agreement. The segments 
(dotted lines) and the error function (dashed lines) fittings are also showed 
in the plot; in the latter, the coincidence with the logistic-like case is so high 
that both curves are almost indistinguishable. 

From each one of the fittings the corresponding expression for 0(t) has been 
estimated. The comparison between them is shown in Figure 5, where we plot 
for simplicity only one of the three cases presented in Figure 4 (the two cases 
non-shown exhibit a very similar behavior). We observe that for the 'error' and 
the 'logistic-like' cases, peaked 0(t) functions with very similar characteristics 
are obtained. The 'segments' case, in turn, leads to a discontinuous expression 
for 0(t) which slightly differs from the other two. So, we can conclude that 
the 'segments' fitting gives a poorer estimate for the behavior of 0(t) and this 
can infiuence the final value for Rq. 

We note that in this specific application for phages a new definition of Rq 



is necessary, as can be seen by inspecting (Q. To this end, we must find the 
equihbrium states of the system 



^ = -PXit)Vit) 

-- = -/3X(t)\/(t)+ / pX(t-t')Vit-t')k{t')(l){t')dt' (20) 

at JO 

and their stabihty. Introducing X{t) = Xgq + 5X{t) and V{t) = \4q + 6V{t) 
and linearizing about the equihbrium states (Xgq, 0) and (0, Veq) one can check 
that the basic reproductive ratio 

Ro = / k{t)(l){t)dt (21) 

Jo 

must be higher than 1 for a successful phage growth. Making use of ( TT^ 

"'-mi w„/' — m— ''" 

which is the burst size. This result simply demonstrates that in the case of 
phage-bacteria interactions the burst size plays the role of a basic reproductive 
ratio (the infection is successful only for Rq > 1). 



7 Conclusions 

In the present paper, we have derived a generalization of the basic model of 
virus dynamics by considering a more accurate life cycle for viruses and cells 
which includes distributed delays for mortality and the eclipse phase. As a re- 
sult, we have showed how the infection dynamics gets modified. As discussed 
above, our main motivation has been to present a rigorous approach to this 
problem, as many times delays have been introduced in this kind of models 
just by intuitive or ad hoc arguments. For this reason, we have provided here a 
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mesoscopic derivation based on explicit balance equations that provide a very 
accurate description of the underlying dynamical process. In our approach, the 
life cycles are implemented in a probabilistic way by the distributions func- 
tions ifxi V^Y, ^v ^"^^ 0- Then, this is a very powerful and general formalism, 
provided that one has the data necessary to evaluate these functions. 

We have carried out a formal analysis of the equilibrium states and their 
stability. Furthermore, we have illustrated how the model works for some 
simple situations of interest. Specifically, for phage-bacteria interactions we 
have been able to provide analytical expressions that may serve to estimate 
the function 0(t) from the one-step growth curves. 

In short, the main conclusions obtained from our study are the following: 

i) The mesoscopic formalism presented here allows to reduce the BMVD of 
3-species to only 2 species (X and V). Then, albeit our model requires a more 
complex mathematical treatment, this simplification can be an interesting 
advantage. 

a) We have formally proved that the stability diagram of the BMVD is in- 
sensitive to any delays considered. It means that the model has always an 
equilibrium infected-free state which becomes unstable for Rq > 1^ which is 
exactly the same condition necessary for the existence of a stable infected 
state. This generalizes similar results previously found (Culshaw and Ruan, 
2000; Nelson and Perelson, 2002; Wang et. al., 2007) that reached the same 
conclusions for more specific cases. 

in) The reproductive ratio Rq and the virus loads are in general very sensi- 
tive to the distributed mortalities considered. It proves one needs to know in 
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detail the cellular life cycles, in special for viruses, to describe the infection 
process. Models based on exponentially distributed death rates will provide 
non-accurate results for the virus dynamics, which can lead to wrong predic- 
tions concerning the success or failure of an immune response or vaccination. 

iv) For phage-bacteria interactions, we have found that fittings of the one-step 
growth based on logistic-like and error functions yield very similar expressions 
for 0(t). From the analysis shown here, it is not possible to determine which 
one of them is more accurate. Anyway, it is clear that both cases give better 
and more realistic estimates for 0(t) and Rq than fittings based on three 
segments. 

In short, we have found that introducing age-distributed processes in the 
BMVD may strongly modify the dynamics of viral infections. These correc- 
tions can be of great interest when the effects of an immune response are also 
considered in the model. Then, the dynamics of the model is expected to be- 
come richer (as happens without delays, too) and the role of the cellular life 
cycles could be more dramatic. Specifically, we expect that age-distributed 
processes can be able to induce new dynamical patterns as periodicity or 
chaos, in the line of recent works on this field (Liu, 1997; Buric et. al., 2001; 
Canabarro et. al., 2004; Wang et. al., 2007). We will address these ideas in a 
forthcoming paper. 



Appendix. Derivation of the model. 

We have introduced ^Pxi't), V'y(it) and i^yit) as the mortality PDF's (see Figure 
1). So that, the probability that a target cell which 'was born' at time t = 

21 



has not died yet at time t is given by $x(^) (hereafter we will refer to it as 
the "survival probability") according to 



^x{t) = l- f Vx{t')^t\ (23) 

and analogous arguments hold for $y(t) and ^^/(t). 
Then, we can write the balance equation for the population densities as 

X(t) =X(0)$x(t)l^x(0, t) + /* Jx{t - t')<^x{t')nx{t - t',t)dt' (24) 

Jo 

y(t)=r(0)$y(t)+ / Jy(t-t')$y(i')dt' (25) 

Jo 

V{t) = V{0)<l>v{t)ftv{0, t) + r Jv{t - t')<l>v{t')Qv{t - t', t)dt', (26) 

Jo 

where Jx(t) represents the density of particles of species X appeared at time 
t, with equivalent definitions for Jrit) and Jv{t). The function Qx(t — t',t') 
is the probability that a particle X do not become infected during the time 
interval (t — t',t), while f2y(t — t',t) is the probability that a virus has not 
been adsorbed by a cell during that interval. So that, the balance equation 
fl2^ simply says that the density of particles X at time t is given by the initial 
density of particles X{0) not infected yet and still alive, plus those target cells 
appeared at any time so far, provided they have neither died nor become 
infected yet. The meaning of Equations 025ti26p follow analogous arguments. 

Regarding the functions Q, their explicit form can be found in the following 

way. For Qx "we take 

rt' 



^x{t — t',t) = exp 



t-t' 



pV{t")df 



(27) 



which corresponds to the solution of the infection equation dX/dt = —(3XV. 
As the infection is considered independent on the other processes (death and 
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production of new cells by the host), the solution of that ODE within the 
interval (t — t', t) gives us a proper definition for the probability Qx{t — t' , t). 
Similarly, from dV/dt = —j3XV we can write 



flvit — t',t) = exp 



px{t")df 



t-t' 



(2J 



The validity of the expressions fl27tl28p can be demonstrated from more rigor- 
ous arguments using the age-structured models by Vlad and Ross (2002). In 
fact, our model (12411261) can be seen as a particular case of the very general 
model by Yadav and Horsthemke (2006), which was in turn based on the orig- 
inal work (Vlad and Ross, 2002). Accordingly, we will follow the formalism in 
(Yadav and Horsthemke, 2006) to derive our model. 

First, we differentiate the system (I24ti26p with respect to t: 



^ = -X(0)nx(0,t) b^(t) +mi)$x(t)] + Jx{t) 

- f Jx{t-t')^^{t')nx{t-t\t)dt' 

Jo 
-l3V{t) f Jx{t - t')^x{t')nx{t - t',t)dt' (29) 

JO 

dY /•* 

— = -Y{0)ipy{t) + Jrit) - Jrit- t')Mt')dt' (30) 

dr JO 

^ = -V{0)Qv{0,t) [ifyit) + PX{t)<l>v{t)] + Mt) 

- I Jvit- t')ipy{t')nv{t - t', t)dt' 
Jo 

-(3X{t) I Jv{t-t')<^v{t')^v{t-t',t)dt' (31) 

Jo 



Then, we introduce ( IMl) and (1261) into ( l29l) and ( 13T]) . respectively, so we obtain 
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dX 

— = -X(0)fix(0, t)vAt) + Jx{t) - (3X{t)V{t) 

'* Jxit - t')^At')nx{t - t', t)dt' (32) 

/o 

dF /■* 

— = -Y{Q)Vy{t) + Jy(t) - / Jy{t- t')ipy{t')dt' (33) 

at JO 

^ = -V(0)fiv(0, t)(^v' W + Jvit) - PX{t)V{t) 

Jv{t - t')'^y{t')VLv{t - t\ t)dt'. (34) 





On the other side, we divide (l24l) by fix(0,t) and transform that equation to 
the Laplace domain (again, we denote the Laplace transform of a function by 
the brackets [■]s with the conjugate variable s). After some easy algebra, it 
can be written as 



1-b 



X\s 



X 



fix(0,t) 



^(0)K]. + K].[^xfixL. (35) 



Finally, introducing the inverse Laplace transform of ( l35l) into (l32l) . the evo- 
lution equation for the species X reads 



dX 



Jx{t)-(3XV~ [ X{t~t')^x{t')nx{t-t',t)dt', (36) 

Jo 



where \l/x is defined in the Laplace domain by ([3]). For the species Y and V 
we can use exactly the same derivation, so that the Equations ( l33|34l) turn 
into 

dV , _ /■* 



dt 
dV 



Jvit) - f Y{t- t')^y(i')dt' (37) 

Jo 

dV /■* 

= jy{t)-(3XV - Vit-t')mv{t')^v{t-t\t)dt' (38) 

dt Jo 

with \&y, \&y defined implicitly in ([3]). 

Hence, we have obtained the general evolution equations (I36tl38p for the model. 
However, note that we still need to give expressions for the densities Jj. From 

24 



Equation ([T]), the number of new target cells appearing at any given time can 
be expressed as 

Jx{t) = A, (39) 

Similarly, the density of infected cells appearing at time t is given by the 
expression 

JY{t)=f3X{t)V{t). (40) 

Finally, the new viruses appeared at time t are given by the function 0(t) (see 
Figure 1) applied over those cell which were infected at any previous time 
t — t', provided they have not died yet. This allows us to write 

Jy(t) = /* Jrit - t')(t){t')<^y{t')dt'. (41) 

Jo 

Once we have the explicit expressions for Jx, Jy and Jy, our model takes the 
final form (12]): 



dX 



clA r 

-— = \-pxv- / x(t-t')^x(t')nx(t-t',t)dt' 

at Jo 

^ = f3XV~ f'Y(t-t')^Yit')dt' 
at Jo 

^ = t (3X(t - t')V(t - t')k(t'U(t')<!>Yit')dt' 
dt Jo 

-I3XV - I V{t-t')^v{t')^v{t-t',t)dt'. 
Jo 
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Figure Captions 

Figure 1. Scheme of the BMVD model with distributed delays for mortality 
and the echpse phase. 

Figure 2. Virus loads obtained numerically from the general model ([2]) for 
the case of Gamma-distributed death times. In the legend we show the values 
of the parameter a = ax = ay = ay used. The other parameters used are 
/3 = 0.02, k = 50, T = 2, T*x = 10/a, t*y = 10/a, r*y = 10/a. 

Figure 3. Definition of the eclipse time, the rise rate and the burst size in 
terms of the one-step growth curve. 

Figure 4. One-step growth for phage T7 inside E. coli. The lines shown rep- 
resent the fittings from the functions in Table 1 (solid lines represent the 
logistic-like fitting, dashed lines correspond to the error function, and the dot- 
ted lines the segments fitting). Symbols correspond to experimental results, 
obtained from You et al. (2002). Circles, triangles and diamonds represent the 
osg-curve for the host growing at 0.7, 1.0 and 1.2 doublings/h, respectively. 

Figure 5. Comparison between the function 0(t) predicted from the three 
different fittings proposed in Table 1. Results shown correspond to the case of 
growing at 0.7 doublings per hour shown in Figure 3. The solid, dashed and 
dotted lines correspond to the predictions from the logistic-like, error function 
and segments, respectively. 
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Table 1 

Characteristics of the three functions proposed for fitting the osg curve, with 
their exphcit expressions for the echpse time, the rise rate and the burst size. 
From (1191) . the estimations for 0(t) are also shown. 



Vit) 



0(t) 



Eclipse Rise Burst 



time rate size 



Three 
segments: 



0; t <r 



H{t-T) 



Vex 



r{t- t); t < t < t + ^ 



V^ 



-H{t-T-^ 



V^: t>T + 



Voc 



Error function: 


Vcx> 

2 


>-erf(^rv/^)' 


r -r^TT{t-Tf /W 
A^ 


T - 


-2/r 


rVoo/4 


Voo 


Logistic-like: 


V^ 


^g-r(t-T) 


T - 


-2/r 


rKo/4 


V^o 


l + e-r(t-r) 


l^l.^r{t-T)Y 



Table 2 



Values of the parameters obtained from the fittings shown in Figure 3 
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Three segments Error function Logistic-like 



Burst size (PFU/ml) 35.9±1.5 (o) 36±1 (o) 37±1 (o) 

35.7±0.9 (A) 37.8±0.6 (A) 38.3±0.6 (A) 

75.2±2.5 (O) 75±25 (O) 76±25 (O) 

Rise rate (PFUml-imin-^) 1.5±0.1 (o) 1.7±0.2 (o) 1.8±0.2 (o) 

3.4±0.4 (A) 3.6±0.2 (A) 3.8±0.2 (A) 

4.9±0.2 (O) 5.9±0.65 (O) 6.3±0.75 (O) 

eclipse time (min) 24.1±0.9 (o) 25.6±0.4 (o) 26.3±1.6 (o) 

21.1±0.5 (A) 21.4±0.1 (A) 21.7±0.4 (A) 

17.9±0.3 (O) 19.1±0.25 (O) 19.5±0.85 (O) 
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